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ABSTRACT 


Starting with certain identities obtained by Reid and Redheffer for 
general matrix Riccati equations we give various algorithms for the case 
of constant coefficients. The algorithms are based on two ideas -- first, 
relate the RE solution with general initial conditions to anchored RE solutions 
and second, when the coefficients are constant the anchored solutions have 
a basic shift-invariance property. These ideas are used to construct an 
integration free superlinearly convergent iterative solution to the algebraic 
RE. Our algorithm, arranged in square-root form, is thought to be numerically 
stable and competitive with other methods of solving the algebraic RE. 


Vi 
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I. INTRODUCTION 


We shall be concerned here with the solution of the following matrix 
Riccatl differential equation (RE) which lies at the heart of many systems 
problems, e.g. optimal control (e.g.[3]) and least-squares estimationCe ,g. [^ ] ) 

— P(t) = FP(t) + P(t)F'^’ + Q - P(t)DP(t), P(0) = Hg 

In this paper F, Q, and D are assumed to be known constant matrices* 

It has been shown in a wide body of literature that if Pg(.) is the 
solution of 

~ Pg(t) = FPg(t) + Pg(t)F'^ + Q - Pg(t)DPg(t), t 2. 3 , Pg(s) = 0 (P) 

i.e, a solution of (1) anchored to zero at time s, then one can recover 

P(.) from P„(.) through the relation 
a 


P(t) = Pg(t) + 4(t,S|Pg(.))[l + P(s)li(t,S|Pg(.))]"''p(3)*^(t,S|P3(.)) (B) 


Where 4> ( . , . | P„( . ) ) and y(.,.|Pg(.)) are defined by 


at 


$(t,3|Pg(.)) = [F - Pg(t)D] «(t,s!P„(.)) , «(S,S|Pg(.)) = I (i^) 


p(t,s|p ( .) ) 




s|P-(.))D*(T,slP_(.))dT (5) 
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Equation (3) has been known at least since the work of Sandor[5] who obtained 
it in a slightly different form (he assumed that P(s) is nonsingular). To 
the best of our knowledge the representation (3), in fact for general nonsym- 
metric REs, is due to Reid, see e.g. [6, eqn. 2.83. Reid's notation for 
the arguments of ‘t' and u has been adopted since it emphasizes the role of 
P_(.). In an estimation context, relations such as (3)-(5) have been discovered 
by Lainiotis[B] via his partitioning approach and also by Worable and Potter[9]- 
[10]. Lainiotis' work highlighted the importance of the smoothing error 
variance relation 


P(s|t) = tP"‘'(s) + y(t,s!P„(. ))]"'' 


( 6 ) 


In this paper we review certain integration-free representations for 
RE solutions and use these recursions to construct a superlinearly convergent 
solution to the algebraic RE. Our results are expressed in algorithmic form. 

II. SHIFT-INVARIANCE PROPERTY OF Pg(.), 4> ( . , . 1 Pg( . ) ) and y(.,|P 3 (.)) 

When F, Q, and D are constant, it has been noted independently by Womble 
and Potter[9]-[ 10] , SidhuC? , section III. 8], and Lainiotis[8] that Pg(-)» 

• j • I ) and u(.,.lPa(-)) have an important shift-invariance property . 

5 o 

namely that for all t ^ s ^ 0, 


Pg(t) = PQ(t - s); ^(t.slPgC.)) = '!>(t-s,o! Pq( . ) ) ; 
u(t,s|P„(.)) = y(t -s,0lPn(.)) 


( 7 ) 


2 
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which prompt the simpler notation 4>Q(t) = 4i(t,0 (Pq( . ) ) and = u(t ,0 |Pq( . ) ) . 

Thus (3) Gan now be written as 

T 

p(t) = pQ(t-s) + $Q(t-s) [I + P(s)ng(t'-s)“^P(s)$Q(t-s) (8) 

an expression that leads to the following recursive algorithm given independently 
by Womble and Potter[9]-[ 10] , Sidhu [7], and Lainiotis[8] : 

Algorithm 1 ; Recursion for P(k5). k = 1.2.3. ... 

(i) For a preselected 6 > 0 Integrate (2), (4)-(5) to obtain Pq( 6), *g(5), uq(6). 

(ii) Compute P(k6), k = 1,2,3i ... recursively, starting with P(0) = IIq, through 

T 

P((k+1)6) = Po(fi) + ^ P(k6)uo(«)]"^P(k5)<tQ(fi) (9) 

Note that step (ii) is an integration- free step. Equation (9) is of 
course just (8) with t = (k+1)6 and s = k6 . The simplicity of (9) makes it an 
attractive candidate for use in numerical computation. However, this algorithm 
has certain undesirable features which include numerical sensitivity to the 
Initial matrix Hq, and a lack of symmetry. One can avoid the influence of IIq 
on the recursions by computing PgCks) and then obtaining P(.) through the 
relation 

P(kfi) = PQ(kd) + $Q(k6)[I + nouo(k6)]"^nQ$o(k6) (10a) 

[this is just (8) with s = 0 and t = ks] or when Hq is nonsingular one could 
use the symmetric relation 

-1 -1 T 

P(k6) = Po(k6) + 4>Q(k5)[nQ + Mo(^'S)] (^0^) 


JPL Technical Memorandum 33-799 


3 


The result (10) is better numerically because it avoids the cumulative effects 
of an ill-conditioned Hq which can propagate through (9). 

Note that (10) requires the computation of <S’Q(k6) and UgOcfi); but these 
terms are of interest. Monitoring Oq gives an on-line measure of stability 
of the time varying linear system (i|); and Pg is useful for analyzing the 
advantages to be gained from smoothing. 

III. RECURSIONS FOR Pg(.), «g(.) AND Pg(.) 

It might seem that the calculation of Pg, 4ig, and Ug, needed in (10), 
would entail considerably more computation than does the recursion (9). 

This is in fact not the case as is shown in the following lemma due to Reid 
[6,p.21] and Redheffer[ 1 1 ] , whose results have here been modified to reflect 
the shift-invariances (7). 

Lemma 1 ; For t ^ s >. 0 we have 

T 

Po(t) = Pg(t - s) + «g(t - s)[I + Pg(s)pQ(t - s)]"^Pg(s)4'g(t - s) , 

( 11 ) 

$g(t) = 4'g(t - S)[I + Pg(s)Pg(t - 3)]"^4’q(s), (12) 

Pg(t) = Pg(s) + $g(s)Mg(t - s)[I + Pg(s)llQ(t - s)]“^«'g(s). (13) 

More recently, this lemma has been studied in independent work by Sidhu 
and Desai[1] and Ljung et al.[2] where it has been shown that Pg(.), *1 ’q(*) 
and Pg(.) can be interpreted as elements of a scattering matrix that can 
be used to solve two-point boundary value problems. Then the results of 
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the lemma have an interpretation in terms of adding layers of a scattering 
medium. The time-varying and nonsymmetrir cases are also considered in 
[l]-t 2 ]. 

By setting t = (k+1)s and s = 6 in these equations one can obtain recur- 
sions that provide Pq(1c 5), ^QCks) and ugCkS) at k = 1 , 2 , 3 ,...; this is left 
to the reader (see also [16]). Instead we present a recursion that is of 
special value when solving the algebraic RE. 

IV. AN INTERVAL DOUBLING ALGORITHM 

An algorithm that recursively generates VqC')* at 

the points t = 2^6, k = 0,1,2, ..., i.e. progresses geometrically (doubling 
the interval at each recursion), follows readily from lemma 1. For this, 
set s = 2*^6 and t = 2s = in the equations of lemma 1. Thus we have: 

Algorithm 2 ; Interval-Doubling. Integration-Free Recursions 

(i) As for algorithm 1, for a preselected $ > 0, compute Pq(5), #0(6), 

yo(a) (2), (^)-(5). 

(ii) Then recursively compute Pq(2^6), «q(2^(S), jiq( 2^6) through the recursions 
Po(2^*-'6) = Pq(2*^6) + $o( 2^6)[I + Pg (2*^6 )pq (2*^6 ) ]“’Po(2*^6 )$J(2‘^6 ) 

(liO 

4o(2^'*'''6) = $q( 2^6)[I + Po(2^6)ug(2^6)]"''$o(2‘^6) (15) 

= Wq( 2^6) + $J(2*^6 )uo(2*"6)[I + Pg (2^6 )uq (2^6 ) j-l $^(2*^6 ) 

(16) 


RFPRODUCmiLir)f OF THE 

.KUUN’AJi PAGE IS POOR 
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When nonsingular matrices, the recursions can be 

put in a symmetric form. The symmetric form involves additional matrix inver- 
sions. However, exploiting symmetry can reduce certain of the computations 
and numerical errors. We omit details of such a formulation because in Section 
V the interval doubling algorithm is discussed using matrix square roots 
which implicitly preserve both symmetry and positivity of ihe variance matrices. 

Algorithm 2 can be used to solve the algebraic RE: 

0 = FP + + Q - PD? (17) 

since P = ^lim P(t). When (F,Q^^^) is stabilizable and (F,D^'^^) is detectable, 
t ® 

the limit is independent of the initial value IIq. Using IIq = 0 it follows that 

P = Pq(®) = , lim Pf,(2*^5) (18) 

The idea of using interval doubling algorithms for solving the steady- 
state problem and hence the algebraic RE was proposed in [12]. There, an 
iterative method, quite different from the one given here, was developed 
for discrete-time Lyapunov and Ricoati equations. Interval doubling is discussed 
in [1]-[2] and, as noted there, doubling methods have been used in radiative 
transfer problems as well . 

It is interesting to note that the above results hold without change 
for the discrete-time Riccati equation, i.e. for i= 0,1,2, ..., 

P(i+1) = 'tP(i)'p'*^ + Q - [4'P(i)H'^+C][l+HP(i)H'^]“‘'[c'^+HP(i)'|j’^] , P(0) = 

(19) 


f 


6 
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In this case, we have 5=1 and the algorithm is initialized with 


PoO) = Q - CC"^; 

V, A SQUARE-ROOT INTERVAL DOUBLING ALGORITHM 

In most applications the RE (2) is such that 5 > 0 can be so chosen 
that Pq(^) is s positive definite matrix. It then follows that Pq0{6) is 
also positive definite for all k > 1; and thus we can arrange the doubling 
algorithm in a symmetric form which is readily put into a square root recursion. 
Thus, using orthogonal transformation techniques (see [I3])i which have been 
successfully applied to least squares estimation [14]-[15], we now develop 
recursions for and where 

Pq ( 2 ^^) = ; UqCbI'S) = ; $0(2^6) = (21) 

Here, Aj^ and Fj^ are by construction upper triangular matrices. The definition 
of Aj^ and Fj^ is of course arbitrary, our asymmetric choice in (21) results 
from the observation that while Po(-) is a variance matrix, infor- 

mation matrix. Further, this choice reduces the number of arithmetic operations 
in our square-root algorithm. 

To make our square-root algorithm transparent , we start by rewriting 

(i^o-de) 


■'Lll'k.! = ■'kl'k + Ck'k>’'tl * Vk3‘’<‘'k*k) (23) 
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where 


The matrix inversions appearing in (22)-(2*l) can be avoided by utilizing 
certain properties of orthogonal transformations: 

Lemma 2 : Let S be an orthogonal matrix chosen so that 


n p n p 

" l[^ i^] S = [Rio] (25) 


where 1^^ is the n x n identity matrix, and S is partitioned consistent with 
(25) as 


REPRODUCIBILITY ()F Tli 
ORIGINAL PAGE Ig POOP 



Then, 


T T 

[I„ + xx^]"'' = SiiS,! ; [Ip + xTx ]“1 = S22S22 ( 26 ) 


Proof I From (25) and the orthogonality of S it follows that Ijj - and 

T -T -T T T 

X = RS 21 • Thus, R = S^i and X = S^iSg-]. Also, from (25) we have [I^ + XX*^] 

= Rr"^, and hence the first of (26) follows immediately. Next, we note that 

■ T T T T 

since S is orthogonal, S 22 S 2 i + Si 2 S-|'^ = 0 and 322^22 ®12®12 " ^p> hence. 
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T X T T -1 ~T 

S22(Ip + X^X)S22 = ^22^22 ^22 ^21^11^11^21 ^22 ” which imnicdi&t'dy 

gives us the other relation in (26). 


Having established this lemma it is easy to obtain the following: 


Algorithm Sauare-Root Interval-Doubling Algorithm 
(i) Construct an orthogonal such that 


[ I„ I r,^,^ ] I 0 ] 


(27) 


where R^^^ is upper triangular. 

(k) ( 

(ii) By using implicitly defined orthogonal transformations and Sj, 

construct upper triangular matrices and such that 


r I T r I o 1 

E • ‘*’kAk^22 ^ ^A ■ ^*k+1 ' ^ 


( k) fT 

where = (S^^ ^ ^k'*'k 
(iii) 


-1 (k) 

‘•’k+l = '*‘k’’k ^11 ^k 






f 




f30; 


Note that significant computational saving can be realized by exploiting 

(k) 

the triangular nature of Aj^, and . Details of the orthogonal trans- 
formation can be found in [13] and [14]. 
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APPENDIX 


COMPUTER MECHANIZATION OUTLINE FOR SQUARE-ROOT INTERVAL DOUBLING 


INPUT : r(N,N) A(N,N) Upper Triangular Matrix Factors 

^(N,N) Transition 

OUTPUT : r, A, (J), updated 


PROGRAM MATRICES: 


A = 


'WORK" 


■ r ‘ 

1 

1 

ie-L. 

1 c 

1 

1 

}» 

_WORK_ 

} 


N 

N 



N 

N 


N 

W = 


i 

1 

1 

1 

!i3' 

”21 

1 ^^22 

1 

1 

1 

1 

1 


MECHANIZATION SKETCH USING QUASI-FORTRAN 


1 

2 


DO 2 I = 1, 2N 
DO 1 J = 1, 3N 
W(I,J)=0. 
W(I,I)=1. 


W = 


DO 3 J = 2N + 1, 3N 
3 W(J - N,J) = 1. 


10 0 


0 10 


W 23 = I 


Tn<*PRODtICIBlLITY OF THE 
oSm. PAOE ta POOR 


4 

5 


DO 5 I = 1,N 
DO 5 J = I,N 
a = 0. 

DO 4 K = I,J 
a = a + r(J,K)* A(N+I,K) 
W(I,J + 2N) = a 


T T 

W^ 3 = = (FA) 

NOTE A(N + I,J) = A(J,I) 


CALL HHL (W, 2N, 2N, 3N, N) 
.T 


I 0 X 
Oil 


HHL 


T 

11 

T 

'21 


^21 

T 

5 

22 


R 
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This is step i, page 9 


6 

7 


DO 7 J = 1,N 
DO 7 I = 1,N 
0 = 0 . 

DO 6 K = 1,1 
0 = O + A(N + I,K)K(j)(J,K) 
W(I,J + 2N) = 0 



8 

9 


DO 9 J = 1,N 
DO 9 I = 1,N 
0 = 0 , 

DO 8 K = 1,N 
0 = 0 + W(N + I, 
A(I,J) = 0 


^ Al - 


N + K>*W(K, 2N + J) 


CALL HHL (A, 2N, N, N, N) 


- 



sL 

a’^ 

L -< 

HHL 

0 




Eq. (28) 


DO 11 I = 1,N 
DO 11 J = I,N 
0 = 0 . 

DO 10 K = I,J 

10 0 = 0 + W(I,K)* r(K,J) 

11 A(I,J) = O 



(upper triangular) 


DO 13 J = 1,N 
DO 13 I = 1,N 
o =0. 

DO 12 K = I,N 

12 0 = O + A(I,K)* <})(K,J) 

13 f(N + I,J) = O 




r (J. 


14 

15 


■ 

DO 15 J = 1,N 
DO 15 I = 1,N 
0 = 0 . 

DO 14 K = 1,1 
0 = 0 + W(K,I)* f(N + K,J) 

A(I,J) = O 


•^1 “ ® 11 ^ 
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CALL UIN(r, 2N, W, 2N, N) 



DO 17 J = 1,N 
DO 17 I = 1,N 
a =0. 

DO 16 K = 1,J 
16 a = a + (j).(I,K)*W(K,J) 
w(i,N+j) = a 





CALL HHU(r, 2N, N, N, N) 


“ - 

r 

HHL 

r 

Y 


0 

— 


L J 


Updated F, Eq . (29) 


DO 19 J = l.N 

DO 19 I = 1,N 

a = 0 . 

DO 18 K = 1,N 

18 a = a + W(K,N+K)* A(K,J) 

19 I) = a 


Updated (j), Eq. (30) 
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SUBROUTINE HHL(W, NRMAX, NR, NC, NOT) 


c ** Lower Triangularization of W * * ( SEE FIGURE ) 
c W(NR, NC) , with NRMAX row dimension (NR .LE. NRMAX) 
c NCT ,LE. NR, NCT ,LE. NC 

c TW computed, with T orthogonal. The result, a partially lower 
c triangular matrix is stored back in the bottom part of W 
c Dimension of low triangular matrix is NCT 


DIMENSION W (NRMAX, NC) 

DOUBLE PRECISION DELTA 
DATA Z/0./, ONE/1./ 

NRPNC = NR + NC @ 'NR' PLUS 'NC' 

NST = NC - NCT + 2 


DO 50 J = NC, NST, -1 
Delta - Z 

JDIAG = NRPNC - J NRMAX 

DO 10 1=1, JDIAG 

10 DELTA = DELTA + W(I,J)**2 
SIG = SQRT ( DELTA ) 

IF (W (JDIAG, J) .GT.Z) SIG = SIG 

FIGURE 

W(JDIAG,J) = W(JDIAG,J) - SIG 
ALFA = ONE/SIG*W(JDIAG,J)) 

IF (J.EQ.l) GO TO 35 
JMl = J - 1 
DO 30 K = 1, JMl 
DELTA = Z 

DO 20 1=1, JDIAG 

20 DELTA = DELTA + W(I,K)*W(I,J) 

DELTA = DELTA*ALFA 
DO 30 1=1, JDIAG 

30 W(1,K) = W(I,K) + DELTA*W(I,J) 


'^ll 

«12 

,^21 

'^22, 


nm 


HHL 


NC 




^21 

V. ' 

”22 


f 


/ 

Unchanged 


35 DO 40 I = 1, JDIAG 
40 W(I,J) = Z 

W (I, JDIAG) = SIG 
50 CONTINUE 
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SUBROUTINE HHU(W, NRMAX, NR, NC, NOT) 


c * * Upper Triangularization of W * * (SEE FIGURE) 
c W(NR,NC) , with NRMAX row dimension (NR .LE. NRMAX) 
c The first NCT columns are triangularized. 
c T W computed, with T orthogonal. The result, a partially 

c upper triangular matrix is stored back in the top part of W. 


DIMENSION W (NRMAX, NC) 

DOUBLE PRECISION DELTA NRMAX 

DATA Z/O./, ONE/1./ 


DO 40 J = 1, NCT 
DELTA = Z 
DO 10 I = J, NR 
10 DELTA DELTA + S(I,J)**2 

SIG = SQRT ( DELTA ) 

IF W(J,J) .GT. Z) SIG = -SIG 
W(J,J) = W(J,J) - SIG 
ALFA = 0NE/(SIG*W(J,J)) 

JPl = J+1 

DO 30 K = JPl, NC 
DELTA = Z 
DO 20 I = J, NR 

20 DELTA = DELTA + W(I ,K) *W(I,J) 

DELTA = DELTA*SIG 

DO 30 I = J, NR 
30 W(I,K) = W(I,K) + DELTA*W(I,J) 

DO 35 I = JPl, NC 
35 W(1,J) = Z 

W(J,J) = SIG 
40 CONTINUE 


NCT NCT 


”11 

^^ 12 ] 

}nr nct / 





HHU ^ ^ 

^21 

^22 



NC Unchanged 

FIGURE 


16 


JPL Technical Memorandum 33-799 



SUBROUTINE UIN(W, NWMAX, WINV, NWINM, N) 


DIMENSION W(NMAX,N), WINV(NWINM, N) 

c W upper triangular 
c WINV = W INVERSE COMPUTED 

c WINV can replace W 

DOUBLE PRECISION a 

c It is good practice to do matrix inversion in double precision 

WINV (1,1) = l./W(l,l) 

DO 20 J = 2,N 
WINV(J,J) = l./W(J,J) 

JM = J-1 

DO 20 K = 1, JM 
a = 0. 

DO 10 I = K, JM 
10 0=0- WINV(K,J*W(I,J) 

20 WINV(K.J) = 0*WINV(J,J) 

RETURN 
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